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Abstract. The weakly interacting trapped Bose gases have been customarily described using the mean-field approximation 
in the form of the Gross-Pitaevskii equation. The mean-field approximation, however, has certain limitations, in particular it 
can not describe correlations between particles. We introduce here an alternative variational approach, based on the correlated 
Gaussian method, which in its simplest form is as fast and simple as the mean-field approximation, but which allows successive 
improvements of the trial wave-function by including correlations between particles. 

INTRODUCTION 

Dilute Bose systems trapped in external fields have been a rapidly developing field since the Bose-Einstein conden- 
sation was observed almost a decade ago. Theoretically the mean-field approach in the form of the Gross-Pitaevskii 
equation 1 1] has been widely and successfully applied to these systems. The computational complexity of the method, 
and thus the computational time, is independent of the number of particles N, in other words it is of the order of 
0(1). Therefore the method can be applied for large (mesoscopic) bosonic systems, especially when combined with a 
pseudo-potential (in the form of the 5-function potential) approximation for the interaction potential between particles. 

However, the mean-field method has certain limitations, in particular it cannot be easily extended to include 
correlations between particles. Correlations become important for systems with higher densities and/or stronger 
effective interactions. Such strong interacting regimes, where the mean-field theory breaks down 1 2], are now routinely 
achieved experimentally by use of Feshbach resonances. 

Rigorous many-body methods, like the diffusion Monte-Carlo method which include all correlations, has 

computational complexity of the order 0(N 2 ) and therefore can only be applied for smaller systems. Again, for 
relatively dilute gases only few simple types of correlation are expected to be important, and including the full 
machinery of rigorous few-body methods is perhaps by far an overkill for these systems. 

Recently, several methods with computational complexity 0(1) have been proposed for finite-range and zero- 
range 1 7] interactions, where the trial wave-function can incorporate two-body correlations. However, these methods 
can not be easily extended to include higher order correlations. 

We introduce here yet another approach which has an important advantage over the existing methods. Namely the 
approach can incorporate any desired number and type of correlations - from an uncorrected wave-function with 
computational complexity of 0(1) at one end, to fully correlated wave-function with computational complexity of 
0(N 2 ) and higher at the other end. Thus, depending upon the problem at hand one has a possibility to negotiate a 
reasonable trade off between the sophistication of the trial wave-function and the computational time. 

For dilute gases only few types of lowest order correlations should be of importance, and it turns out that for these 
types of correlations the method is yet of 0(1) order of complexity. 

The method is based on correlated Gaussians and amount to a judicious choice of the nonlinear parameters of the 
basis Gaussians. 

METHOD 
Jacobi coordinates 

Consider a system of N particles with masses mi, coordinates r,, i = 1..N, and the Hamiltonian 

N h 2 d 2 



where Vij is the potential between particles i and j and V ext is the external confining potential (a trap). Usually the trap 
is assumed to be harmonic, 

N j 

V ext = E ^a 2 rf • (2) 



It is of advantage to introduce the scaled coordinates, q, = J -^-r,, where m is an arbitrary mass scale. Indeed the 
kinetic energy operator T and also the harmonic trap potential V ext have then a more symmetric form, 

The Jacobian of the transformation from r to the scaled coordinates q is equal to 

f^=n(-) 3/2 - w 

d{r\..r N ) Y\m/ 

If all particles have the same mass m, there is no difference between coordinates r and q. 
One can make a further suitable linear transformation to a new set of coordinates, 

*i = Y,Uij<lj, (5) 

j 

or, in matrix notation x = f/q, where the matrix U is independent of q. The new system of coordinates is called Jacobi if 

i) one of the coordinates, say the Mh, is proportional to the center of mass coordinate R of the system, = y^^R; 

ii) the other N — I coordinates are translation invariant; and iii) the transformation preserves the "diagonal" form 
of the kinetic energy operator. 

The last property implies that the transformation (0 and also any transformation between different Jacobi coordi- 
nates is unitary, UU T = 1 (where T denotes transposition), with the corresponding Jacobian being equal to unity. The 
unitarity means that the so-called hyper-radius p, defined as p 2 = £q?, is invariant under these transformations, 

p 2 -M = K = ^E-^. (6 ) 

With Jacobi coordinates the center of mass coordinate decouples and the hyper-radius it therefore often defined 
without the contribution from the center of mass coordinate x^, 

P 2 = E x < =£q?-iVR 2 . (7) 

i<N i 

One of the possible choices of the Jacobi coordinates is 

x/=iJv = \/^(R/-r i+ i) , (8) 
V m 

where R,- is the coordinate of the center of mass of the first i particles, = 0, and /I, is the reduced mass 

# = ITT ' (9) 

where M, = £^ =1 m k . 

In the following we shall only consider identical particles with m, = m. 

Hyper-radial approximation 

Non-interacting bosons in a harmonic trap 



Let us consider a system of non-interacting bosons in a harmonic trap. This should be a good first approximation to 
a system of weakly interacting bosons in a trap which is smooth at the bottom and spherically symmetric. 



The ground-state wave-function *P of a system of non-interacting bosons is a product 

^=rivofe), (io) 

i 

where V / b(?) is the lowest (s-wave) single-particle state of the trap. If the trap is harmonic, V / o(?) is a Gaussian, 

12 —1/2 

Vb(?) 0,1 e~T- a ° q ~ , where a is the (scaled) oscillator length, and the ground-state wave-functions simplifies to a 
single Gaussian depending only on the hyper-radius p, 

*i> = Y[yo{ qi ) ~Y[e-^ =e-? a °Zi<it =e -i«oP 2 . (11) 

i i 

1 2 

A single Gaussian e~2 a oP is thus an exact solution for a system of non-interacting bosons in a harmonic trap. 
Generally speaking a function of hyper-radius will provide an exact solution to the many -body system in cases where 
the potential energy of the system depends only on the hyper-radius. The harmonic trap is precisely this type of 
potential. 



Weakly interacting bosons 

If the particles in the trap interact only weakly one can assume, following the ideas from the mean-field theory, 
that the inter-particle interactions will effectively lead to a certain modification of the field. The solution will then be 
some square-integrable function of hyper-radius, <$>hr{p), which can be represented as a linear combination of, say, n 
Gaussians, 

*hr{p) - t C s e-^P 2 = YJ[C s e--- a «l (12) 

,V= 1 S ( 

where C s are variational parameters, and the range parameters a s (s — 1. .n) are assumed to be fixed and chosen to 
span the necessary functional space. This trial wave-function is called a hyper-radial approximation. In practice the 
parameters a s are chosen and then optimized in a stochastic procedure using the ideas from the stochastic variational 
method |9]. 



Hyper-radial vs. mean-field 

The variational mean-field approach is based on an assumption that a product wave-function can provide a good 
description of an interacting system. The trial wave-function ^mf is taken as a product of single-particle functions l//, 

^mf =l\v(qi), (13) 

i 

where the functional form of \]/(q) is varied to reach the minimum of the expectation value of the Hamiltonian. 
Assuming that y/ is a square integrable function, one can represent it as a linear combination of Gaussians, 

^)=£c,e->^, (14) 

s 

where the coefficients c s are the variational parameters. The trial mean-field wave-function then becomes 

v aff =n£ c ,«-w j d5) 

i s 

which should be compared with the hyper-radial trial wave-function 

^HR{p)=YX\C s e-^-. (16) 

s i 

The two trial functions (II 51 and d!6i are similar but not equivalent since the sum and the product operators generally 
do not commute. Note that the hyper-radial variational parameters C s are linear, while the mean-field parameters c s are 



non-linear 1 . In practice, however, as we shall show by numerical calculations, both trial functions give rather similar 
results. 

Both functions are totally symmetric and thus do not require an explicit symmetrization. The computational time 
for the variational minimization of the Hamiltonian with both functions is independent of the number of particles. 

The hyper-radial function has an advantage that the center of mass motion can be easily decoupled by a (unitary) 
transformation to relative Jacobi coordinates. Again, the mean-field function cannot be easily improved, while the 
hyper-radial function is only the basis for further improvements. 



Correlations 

Two-body correlations 

The correlation between a pair of particles can be described by a basis function in the form 

<p n=e -^ a P 2 -^^-^\ (17) 



where there are now two independent parameters, a and /3 . The trial wave-function is then a linear combination of 
<E>i2's with different parameters a and j3, 

«P = Y^C^e-^P 2 -^"^-^ 2 , (18) 

.v.;/ 

where C.„, are linear variational parameters. The nonlinear parameters a and j3 are again chosen and optimized 
stochastically. 

The basis function is no longer automatically symmetric over all permutations. It has to be symmetrized with respect 
to particles number 1 and 2 and therefore the symmetrization operator, S, has to be included when calculating matrix 
elements, 

N ^ 1 



S<t>n=[ j ) E*«- < 19 > 



ij 



This is the same type of Faddeev-like decomposition of the wave-function as used in §01. 

Fortunately, only a finite number of different terms appear in calculations of the matrix elements, and the compu- 
tational time is therefore still independent of the number of particles. Indeed the kinetic energy and the external field 
operators are fully symmetric and therefore the explicit symmetrization of the wave-function is not needed for their 
matrix elements. The matrix element for the inter-particle potentials reduces to a finite number of terms, 

N (<t>n\Li<jV u S\<t> l2 )= (20) 



(*12 I (Vi2 + 2(jV-2)V 13 + (W ' 2) f' 3) V 34 ) | 4> 12 ) 
+2(iV-2)<<Di 2 | 
(V12 + V l3 + V 2i + (N- 3)(Vi4 + V 24 + V 34 ) + (jV ' 3) 2 (Af ' 4 V 45 ) 

1*13} 

+ p£a-i-2(tf-2))<* 12 | 

(V12 + 4V 13 + V 24 + V34 + 2 (N - 4) (V15 + V35) + {N - A) 2 {N - 5) V 56 

1*34) 

Each individual matrix element in this expression is readily calculated using the expression J28i in the appendix. The 
structure of the expression basically corresponds to that of |8] where hyper- spherical coordinates were used instead 



indeed the Gross-Pitaevskii mean-field equation is non-linear. 



of the Jacobi coordinates used here. Hyper- spherical coordinates allow an easy implementation of a powerful hyper- 
spheric adiabatic expansion method but, on the other hand, do not allow an easy implementation of higher order 
correlations. 



Three-body correlations 
The three-body correlations can be accounted for by a basis function of the form 

Phi n?> = e-s^-sPfoi-*) 2 -!**-*) 2 , (21) 

where a, j3 and 7 are independent parameters. The trial wave-function is then a linear combination of Om's with 
different parameters a and j3 and 7, 

vp = Q, Hiie -i«iP 2 -iA,(q 1 -q 2 ) 2 -^rv(qi-q3) 2 j (22) 

s,u,v 

where C suv are linear variational parameters, and where the nonlinear parameters a, j3 and 7 are again chosen and 
optimized stochastically. 

This function must be explicitly symmetrized with respect to particles 1, 2, and 3. This symmetrization again results 
in a finite number of different terms as it did for two-body correlations. There are in total 34 different terms and it is 
therefore not practical to write them down here. The computer program can easily catch the identical terms and thus 
reduce the computational complexity down to the order of 0(1), that is, independent of the number of particles. 



NUMERICAL ILLUSTRATIONS 
The Bose system 

We use 87 Rb condensate parameters corresponding to fixed scattering length a s = 100 a.u. and trapping frequency 
(0 = 2n x 77.87 Hz, and vary the number of atoms N = 10 1 — 10 4 . In all cases, the inverse square root of the nonlinear 
parameters fa and 7- are optimized from the random value interval [10~ 4 £>,; lQb t ] (where b t — y/h/{mco) ~ 23095 a.u. 
is the trap length), while for the parameters a the interval was [b t ; I0 i b t ]. In practice only one parameter 0Cq was needed 
to achieve the chosen accuracy goal of three digits on the interaction energy per particle. 

The mean-field validity condition, na\ <C 1, where n is the particle density, is fulfilled for all values of N. Therefore 
the Gross-Pitaevskii results from the literature should be quite accurate and we shall use them as the reference point. 
The other regime, na\ ^> 1, shall be investigated separately. 



Two-body potentials 

We consider only dilute bosonic systems where the properties largely depend upon the low-energy/large-distance 
properties of the two-body interaction, that is the s-wave scattering length a s . In this regime a zero-range pseudo- 
potential given by a delta function, 

V s (r) = 4 -^S(r), (23) 
m 

is proven to provide within a mean-field theory a good approximation to the energy of the system. Applying the 
delta-function interaction with a Hilbert space of a beyond-mean-field theory, however, requires an appropriate 
renormalization |7]. The physical scattering length in d23i should be substituted by its first-order Born approximation 
of the given finite-range potential. 

We shall use the delta-function potential for calculation with the uncorrelated hyper-radial trial wave-function. 

For correlated calculations we shall use four different finite-range potentials of the form 



V(r) =V Q e- rl l h2 + U (i e- rl l <2 , 



(24) 



TABLE 1. The parameters (in atomic units) of the finite-range Gaussian 
two-body potentials of the form V(r) = Vo e ~ r '* + Uge ' c use d in the 
calculations. Nt, is the number of bound states in the potential. The s-wave 
scattering length a s is equal 100 a.u. for all potentials. 



Designation 


b 


Vo 


c 


U 


N h 


H (hard) 


58.69 


1.906 x 10" 7 











S (soft) 


550.0 


1 x 10" 11 











A (attractive) 


10 


-1.906 x 10- 7 








1 


W (well) 


4.4 


5.566 x 10" 5 


10 


-1.125 x 10" 6 


1 



TABLE 2. The interaction energy per particle, § — |/ift), where E is the total energy, 
for the system described in the text. Results are given for the hard-core (H) and soft- 
core (S) potential from TableQwith different trial wave functions (lb - uncorrelated, 
2b - two-body correlations, 3b - three-body correlations) as well as for the 8 -function 
potential with uncorrelated wave-function. The last column shows the Gross-Pitaevskii 
(mean-field) results from 1 10] and | IX]. 





hard-core potential 


soft- 


core potential 


5-function 




N 


lb 


2b 


3b 


lb 


2b 


3b 


lb 


GP 


10 


.329 


.0155 


.0154 


.0179 


.0154 


.0154 


.0154 


.0154 


20 


.599 


.0326 


.0325 


.0373 


.0320 


.0320 


.0320 


.0320 


50 


1.18 


.0832 


.0828 


.0923 


.0795 


.0794 


.0798 


.0792 


100 


1.83 


.165 


.164 


.177 


.153 


.153 


.153 


.151 


1000 


6.29 


1.32 


1.32 


1.09 


1.00 


.999 


.978 


.930 


5000 


13.2 


4.48 


4.47 


2.88 


2.75 


2.75 


2.64 


2.45 


10000 


17.8 


7.27 


7.26 


4.15 


4.02 


4.02 


3.83 


3.58 



where the parameters of the potentials are specified in Tabled The first potential, marked H, is a hard repulsive core, 
the second, S, is a soft repulsive core, the third, A, is an attractive well, and the fourth, W, is a semi-realistic well with 
a repulsive core and an attractive pocket. All potentials have the same scattering length, a s = 100 a.u., and in the dilute 
regime should therefore provide identical energies if correlations are appropriately included. 



Results 

The results are collected in Tables|2]and[3] where we show the interaction energy per particle, § — jhco (where E is 
the total energy of the system), for different combinations of numbers of particles, potentials, and trial wave-functions. 
The absence of a number for the attractive and realistic potential means that there are many strongly bound (collapsed) 
states and an analog of the condensate state located in the trap does not exist. 

TABLE 3. The same as Table|2|for the attractive (A), and realistic 
(W) potentials from Table Q For larger number of particles and 
higher correlations the potentials produce a large number of strongly 
bound (collapsed) states and thus no condensate state could have 
been traced. 





attractive potential 


realistic potential 




N 


lb 


2b 


lb 


2b 


3b 


GP 


10 


-.0021 


.0147 


.0383 


.0154 


.0154 


.0154 


20 


-.0044 


.0264 


.0599 


.0320 


.0320 


.0320 


50 


-.0114 


.0228 


.188 


.0804 


.0802 


.0792 


100 


-.0233 


-.0042 


.344 


.156 


.155 


.151 


1000 






1.78 


1.07 




.930 


5000 






4.33 


3.27 




2.45 


10000 






6.11 


5.09 




3.58 
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FIGURE 1. Interaction energy per particle as function of the number of particles N for the uncorrelated trial wave-function. 



The results for different potentials with the uncorrelated hyper-radial trial wave-function are given in Tables l2l3l and 
also represented on FigQ 

Importantly, the combination of delta-function pseudo-potential with the uncorrelated hyper-radial wave-function 
give results within a few per cent of the mean-field theory. The pseudo-potential therefore seems to be equally well 
suited for both mean-field and hyper-radial approximations. 

One can use this very fast uncorrelated pseudo-potential approximation to a great effect as a tool to optimize the 
parameters of the Gaussians to be used in the more demanding correlated calculations with finite-range potentials. 

The finite-range potentials show large deviations since the uncorrelated wave-function is not suited for them. The 
hard-core potential, as could be expected, is especially bad for the uncorrelated wave-function. The attractive potential 
produces for larger number of particles a strongly bound (collapsed) ground-state and is therefore not shown on the 
figure. 



The results with the two-body correlated trial wave-function are given in Tables l2l3l and also represented on Fig. [2] 
Apparently, inclusion of two-body correlations dramatically improves the results. This seems to support the assertion 
in that the two-body correlations are of the utmost importance for the dilute gases. 

The hard-core potential, although doing much better with the two-body correlated wave-function, is still the farthest 
off especially for large number of particles. The soft-core potential on the other hand is now very close to the mean- 
field results. 



We do not show a separate figure for the three-body correlations as they turn out not to produce large effects on the 
energies apart from potentials with attraction, where the three-body correlations quite expectedly straight away lead to 
a large number of strongly bound (collapsed) states. 

Thus, for model repulsive finite-range potentials and dilute systems the three-body correlations are of much less 
importance that two-body correlations. 
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FIGURE 2. Interaction energy per particle as function of the number of particles A' for the trial wave-function with two-body 
correlations. 



We have introduced a new approach, based on correlated Gaussian method, to investigate dilute Bose systems. The 
approach allows to include consecutively correlations of different orders in the trial wave-function. In its lowest 
(uncorrelated) order with zero-range pseudo-potentials the approach is comparable to the mean-field (Gross-Pitaevskii) 
theory. 

We have performed an exploratory numerical investigation of two- and three-body correlation in a dilute Bose 
system with different number of particles and different finite-range potentials. For the condensate state the two-body 
correlations turn out to be by far the most important and suffice to provide a quantitative description of the system with 
soft-core potentials. 



H. H. S0rensen would like to thank Christoffer Dam Bruun for many interesting discussions and valuable assistance 
during implementation of the method. 



The trial wave-function is represented as a linear combination of correlated Gaussians, |A), which have the form 



where A is a positively definite symmetric matrix and x is a set of (scaled Jacobi) coordinates. Correlated Gaussians 
form a full basis since any square-integrable function can be represented as a linear combination of Gaussians with 
arbitrary precision. The elements of the parameter matrix A can be optimize using the stochastic method |9|. 
The important matrix elements which are used in the calculations are the overlap of two Gaussians, 



CONCLUSION 
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APPENDIX: CORRELATED GAUSSIAN METHOD 




(25) 




(26) 



the matrix element of the kinetic energy operator, 



h 2 



M-toL^M^MlA+tT^W*), (27) 

and the matrix element of the two-body potential V (r; — r 7 ), 

d 3 rV (r)(A | 8(bfjx - r) | A 1 ) = G C; . [V] (A | A'), (28) 

-oo 

where r, — r 7 = bf-x 7 cjj 1 — bfj(A +A f )~ l bij, and G C [V] is the Gaussian transform of the potential 

G C [V] = (^) 3/2 / d 3 rV(r)e-^ r \ (29) 

Other useful integrals 

(A | x T Bx | A') = 3tr((A +A')~ 1 fi) (A | A') ; (30) 
{A\8(b T x-q)\A')= (J^j ' e~^' 2 {A \ A') , where j3 _1 = b r (A +A') _1 b ; (31) 

G t [i]=2^; (32) 

3/2 



[*W1 = (^) ; (33) 



t^i _ ( c \ 



3/2 



Gc[S ~ 1 = 1 ^ 1 • (34) 
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